function gradient = quat2expmapGradient(expmap)

gradient  = zeros(4, 3);

minAngle = 1.0e-7;
v = expmap;
theta = sqrt(v(:)' * v(:));
cosp = cos(0.5 * theta);
sinp = sin(0.5 * theta);

for i = 1:3
   if theta < minAngle
       Tsinc = 0.5 - (theta * theta / 48.0);
       vTerm = v(i) * ((theta * theta / 40.0) - 1.0) / 24.0;

       gradient(1, i) = v(1) * vTerm;
       gradient(2, i) = v(2) * vTerm;
       gradient(3, i) = v(3) * vTerm;
       gradient(4, i) = -0.5 * v(i) * Tsinc;

       gradient(i + 1, i) = gradient(i + 1, i) + Tsinc;
   else
       ang = 1.0 / theta;
       sang = sinp * ang;
       ang2 = ang * ang * v(i);
       cterm = ang2 * (0.5 * cosp - sang);

       gradient(1, i) = cterm * v(1);
       gradient(2, i) = cterm * v(2);
       gradient(3, i) = cterm * v(3);
       gradient(4, i) = -0.5 * v(i) * sang;

       gradient(i + 1, i) = gradient(i + 1, i) + sang;
   end
end
